library(tidyverse)
library(cluster)
library(sandwich)
library(lmtest)
library(MASS)
library(stargazer)
library(xtable)
library(fixest)
library(reshape2)
library(broom)
library(ggridges)
library(RColorBrewer)
library(nnet)
library(margins)
library(sf)

options(error = NULL)

dat <- read_csv("Mindanao_survey_results_Koyama.csv") %>%
  mutate(income_2 = as.numeric(income),
         peaceful = as.numeric(reb_violence_2019 == "none"),
         immigrant = as.numeric(pob != "same"),
         male = as.numeric(gender == "male"),
         othermun = as.numeric(pob == "othermun"),
         otherprov = as.numeric(pob == "otherprov"),
         noBARMM = as.numeric(pob == "noBARMM"),
         maguindanao = as.numeric(language == "Maguindanao"),
         tagalog = as.numeric(language == "Tagalog"),
         iranun = as.numeric(language == "Iranun"),
         urban = as.numeric(Locale == "Urban"),
         president = as.numeric(president_code == 1),
         pm = as.numeric(pm_code == 1),
         capital = as.numeric(capital_code == 1),
         overelement = as.numeric(education %in% c("comp_college", "comp_elementary", "comp_high", "comp_vocational",
                                                   "post_college", "some_college", "some_high", "some_vocational")),
         overhighschool = as.numeric(education %in% c("comp_college", "comp_high", "comp_vocational",
                                                "post_college", "some_college", "some_vocational")),
         witness_rido = ifelse(witness_rido == "Yes", 1, ifelse(witness_rido == "No", 0, NA)),
         witness_reb = ifelse(witness_reb == "Yes", 1, ifelse(witness_reb == "No", 0, NA))) %>%
  mutate(richer_median = as.numeric(income_2 > median(income_2, na.rm = T)),
         immigrant_5 = immigrant * (years_live < 5),
         immigrant_10 = immigrant * (years_live < 10 & years_live >= 5),
         immigrant_long = immigrant * (years_live > 9),
         knowledge = president + pm + capital)

## Support for peace agreement has an NA code "97". Replace them with NA
dat$m1q4[dat$m1q4 == 97] <- NA
dat$m1q5[dat$m1q5 == 97] <- NA
dat$m1q6[dat$m1q6 == 97] <- NA

## A categorical voting variable for plotting purpose
dat$vote_cat <- ifelse(dat$vote == 0, "NonVoter", ifelse(dat$UBJP_vote == 1, "UBJP", "NotUBJP"))

##### Mapping sampled barangays

## Load the white map for Maguindanao
map <- readRDS("gadm36_PHL_3_sp.rds")
map <- subset(map, map$NAME_1 == "Maguindanao")

## Combined strings of municipality and barangay names are used as ID for merger
dat$mun_bar <- paste(dat$Municipality, dat$Barangay)
location <- unique(dat[,c("mun_bar", "Locale")])

## South Upi is referred to by Upi in the map data. Adjust it.
map$NAME_2[map$NAME_2 == "Upi"] <- "South Upi"
map$mun_bar <- paste(map$NAME_2, map$NAME_3)

map$sample <-  map$mun_bar %in% location$mun_bar
map@data <- merge(map@data, location, by = "mun_bar", all.x = T)
map@data$Locale[is.na(map@data$Locale)] <- "NotSample"

# Define Colors and Breaks for Categories
cols <- c(0,4,3)
names(cols) <- c("NotSample", "Urban", "Rural")
# Match the colors to the data
map$col <- cols[map$Locale]

rownames(map@data) <- as.character(seq_len(nrow(map@data)))

plotmap <- st_as_sf(map)

## Figure 1
plot(plotmap["Locale"], col = plotmap$col, main = "Sampled Barangays")

### Evaluation of BTA

## Subset the data to the control group only
## (Please see the appendix for the experimental component of the survey module)
stability <- dat$m3q1[dat$m3g == "C"]
localpart <- dat$m3q2[dat$m3g == "C"]
economic <- dat$m3q3[dat$m3g == "C"]
corruption <- dat$m3q4[dat$m3g == "C"]
postpone <- dat$m2q[dat$m2g == "C"]

btaeval <- data.frame(question = rep(c("Stability", "Participation", "Economy", "Anti-Corruption", "BTA_Extension"), each = 175),
                      response = c(stability, localpart, economic, corruption, postpone))
btaeval$likert <- factor(btaeval$response, labels = c("VeryHigh", "High", "Middle", "Low", "VeryLow"))

custom_colors <- c("VeryHigh" = "#018571",
                   "High" = "#80cdc1",
                   "Middle" = "#f5f5f5",
                   "Low" = "#ff7f7f",
                   "VeryLow" = "#d73027")
btaeval$question <- factor(btaeval$question, levels = c("Stability", "Participation", "Economy", "Anti-Corruption", "BTA_Extension"))

## Figure 2
ggplot(btaeval, aes(x = question, fill = likert)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = custom_colors) +
  labs(x = "Question", y = "Percentage", fill = "Response") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
## To save the graph, use the following
# ggsave("BTAeval_likert.jpg", width = 6, height = 4)

## Support for peace clauses

tjrc <- dat$m1q4[dat$m1q4g == "C"]
shariah <- dat$m1q5[dat$m1q5g == "C"]
decommission <- dat$m1q5[dat$m1q6g == "C"]
paeval <- data.frame(question = rep(c("TJRC", "IslamicCourt", "Decommission"), each = 140),
                     response = c(tjrc, shariah, decommission))
paeval$likert <- factor(paeval$response, labels = c("VeryHigh", "High", "Middle", "Low", "VeryLow"))

## Figure 3
ggplot(paeval, aes(x = question, fill = likert)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = custom_colors) +
  labs(x = "Question", y = "Percentage", fill = "Response") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
# ggsave("policyeval_likert.jpg", width = 6, height = 4)

## satisfaction index
dat$satisfy <- dat$m3q1 + dat$m3q2 + dat$m3q3 + dat$m3q4
dat$satisfy_index <- -scale(dat$satisfy)

## Distribution of priority issues
issue <- as.data.frame(table(dat$pol_issue))
names(issue) <- c("issue", "count")

issue <- issue %>%
  arrange(desc(count)) %>%
  mutate(prop = count / sum(count) * 100) %>%
  mutate(ypos = cumsum(prop) - 0.5 * prop)
issue$issue <- str_to_upper(issue$issue)

issue$ypos <- c(200, 480, 595, 640, 670, 700)
#" Create the pie chart
library(ggrepel)
ggplot(issue, aes(x = "", y = count, fill = reorder(issue, count))) +
  geom_bar(stat = "identity", width = 1, color = "white") +
  coord_polar("y", start = 0) +
  theme_void() +
  geom_label_repel(aes(y = ypos, label = paste0(round(prop, 1), "%")), color = "white", size = 5) +
  labs(fill = "Issue") +
  theme(legend.position = "right") +
  scale_fill_brewer(palette = "Set2") +
  guides(fill = guide_legend(reverse = T))
# ggsave("issue_pie.pdf", width = 6, height = 4)

## The following code creates an interactive pie chart, which is used in the paper.
## The results are the same as above, just different appearance.
## If you want to replicate it, replace "0" in the next line with "1" and run the code.
if(0){

library(plotly)
library(htmlwidgets)
library(webshot2)

issue <- issue %>%
  arrange(desc(count)) %>%
  mutate(prop = count / sum(count) * 100)

# Customize the order of issues
custom_order <- c("PEACE", "ECONOMY", "EDUCATION", "HEALTH", "DISASTER", "LAND")

# Reorder the issue dataframe according to the custom order
issue$issue <- factor(issue$issue, levels = custom_order)
issue <- issue %>% arrange(match(issue, custom_order))

custom_colors <- c("#DDA0DD", "#87CEEB", "#98FB98", "#FFD700", "#E9967A", "#FFA07A")

# Create interactive pie chart
fig <- plot_ly(issue, labels = ~issue, values = ~count, type = 'pie', textinfo = 'label+percent',
               insidetextorientation = 'radial', hoverinfo = 'label+percent+value',
               marker = list(colors = custom_colors)) %>%
  layout(title = "",
         showlegend = FALSE,
         margin = list(l = 50, r = 50, b = 50, t = 50),
         font = list(size = 25))

htmlwidgets::saveWidget(fig, "temp_plot.html", selfcontained = T)
webshot2::webshot("temp_plot.html", file = "issue_pie_chart.jpg", cliprect = "viewport")
}

## Create binary variables for different job categories
expand_categorical_to_binary_keep_na <- function(data, categorical_var) {
  data_copy <- data
  # Create a new column to indicate NA values in the categorical variable
  data_copy[[paste0(categorical_var, "_is_na")]] <- is.na(data_copy[[categorical_var]])
  data_copy[[categorical_var]] <- as.character(data_copy[[categorical_var]])
  data_copy[[categorical_var]][is.na(data_copy[[categorical_var]])] <- "NA"
  # Create a model matrix for the categorical variable
  model_matrix <- model.matrix(~ . - 1, data = data_copy[, categorical_var, drop = FALSE])
  
  binary_vars <- as.data.frame(model_matrix)
  colnames(binary_vars) <- gsub(paste0(categorical_var, ""), "", colnames(binary_vars))
  new_data <- cbind(data, binary_vars, data_copy[[paste0(categorical_var, "_is_na")]])
  
  return(new_data)
}
dat <- expand_categorical_to_binary_keep_na(dat, "job")

## Descriptive summary of covariates

ratio_table <- function(x){
  count <- table(x)
  count <- c(count, sum(count))
  ratio <- count / (sum(count)/2)
  out <- data.frame(t(rbind(count, ratio)))
  out <- out[order(out$count, decreasing = T),]
  out <- out[c(2:nrow(out), 1), ]
  return(out)
}

## Tables 1 and 2 in the Online Appendix are built from the following results

ratio_table(dat$gender)
mean(dat$income_pc, na.rm = T)
sd(dat$income_pc, na.rm = T)
mean(dat$age, na.rm = T)
sd(dat$age, na.rm = T)
ratio_table(dat$job)
ratio_table(dat$education)
ratio_table(dat$language)
ratio_table(dat$Locale)
ratio_table(dat$pob)
ratio_table(dat$reb_violence_2019)
ratio_table(dat$reb_violence_2014)
ratio_table(dat$witness_reb)
ratio_table(dat$witness_rido)
ratio_table(dat$displace)
ratio_table(dat$clan_affiliation)
ratio_table(dat$vote_cat)
ratio_table(dat$president)
ratio_table(dat$capital)
ratio_table(dat$pm)

## Some additionjob_other## Some additional variables for the analysis
dat$cotcity <- as.numeric(dat$Municipality == "Cotabato City")
dat$income_2

## Household income per capita
dat$income_pc <- dat$income_2 / dat$household_size
## Binary indicator if the respondent is richer than the median in terms of the household income per capita
dat$richer_median_pc <- as.numeric(dat$income_pc > median(dat$income_pc, na.rm = T))

## Older than the median
dat$old <- as.numeric(dat$age > median(dat$age))
## Have better knowledge?
dat$better_knowledge <- as.numeric(dat$knowledge > 1)
## Indicator if voted but not for UBJP
dat$vote_other <- as.numeric(dat$vote - dat$UBJP_vote)

## Regression analysis for Figure 5

indx1_full <- lm(satisfy_index ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                   overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                   displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)

## Function for coefficient plot
coefplot <- function(fit, exclude = NULL, term = NULL, order = NULL){
  tid <- tidy(fit, conf.int = T)
  tid_p<- subset(tid, !tid$term %in% exclude)
  if(!is.null(term)){
    tid_p$term <- term
  }
  tid_p$col <- "#555555"
  tid_p$col[tid_p$conf.low > 0] <- "#FF7F50"
  tid_p$col[tid_p$conf.high < 0] <- "#3498DB"
  
  if(is.null(order)){
    ggplot(tid_p, aes(x = estimate, y = reorder(term, estimate))) +
      geom_point(aes(color = col), size = 3) +
      geom_errorbarh(aes(xmin = conf.low, xmax = conf.high, color = col), height = 0.2) +
      geom_vline(xintercept = 0, linetype = "dashed") +
      theme_minimal() +
      labs(title = "",
           x = "Coefficient",
           y = "Variable") +
      scale_color_identity()
  }else{
    tid_p$order <- order
    ggplot(tid_p, aes(x = estimate, y = reorder(term, order))) +
      geom_point(aes(color = col), size = 3) +
      geom_errorbarh(aes(xmin = conf.low, xmax = conf.high, color = col), height = 0.2) +
      geom_vline(xintercept = 0, linetype = "dashed") +
      theme_minimal() +
      labs(title = "",
           x = "Coefficient",
           y = "Variable") +
      scale_color_identity()    
  }
}

exc_m3 <- c("m3gT1", "m3gT2", "m3gT3", "(Intercept)")
fullterm <- c("Male", "Richer", "Older", "Homemaker", "Selfemploy", "Farmer", "Worker", "ElementarySchool", "Highschool",
              "Maguindanao", "Tagalog", "Iranun", "Urban", "Cotabato City", "Immigrant", "Peaceful",
              "Witness Insurgency", "Witness Rido", "Displacement", "Clan Affiliation", "Voted",
              "Voted for UBJP", "Political Knowledge")
## Manual reordering of the graph
fullord <- c(6, 9, 7, 12, 16, 19, 8, 10, 13, 22, 11, 5, 14, 18, 2, 17, 20, 3, 15, 4, 1, 23, 21)

## Figure 5
coefplot(indx1_full, exc_m3, fullterm)
# ggsave("coefplot_full.jpg", width = 6, height = 4)

# Summarize the data to calculate mean and confidence intervals
vote_cat_sum <- dat %>%
  group_by(vote_cat) %>%
  summarize(
    mean_satisfy_index = mean(satisfy_index, na.rm = TRUE),
    ci_low = mean(satisfy_index, na.rm = TRUE) - qt(1 - (0.05 / 2), df=n()-1) * sd(satisfy_index, na.rm = TRUE) / sqrt(n()),
    ci_high = mean(satisfy_index, na.rm = TRUE) + qt(1 - (0.05 / 2), df=n()-1) * sd(satisfy_index, na.rm = TRUE) / sqrt(n())
  )

mean_data <- dat %>%
  group_by(vote_cat) %>%
  summarize(mean_satisfy_index = mean(satisfy_index, na.rm = TRUE))

## Create the ridge plot with points and vertical lines for means, aligning colors
## Figure 6
p <- ggplot(dat, aes(x = satisfy_index, y = vote_cat, fill = vote_cat)) +
  geom_density_ridges(scale = 0.9, alpha = 0.5) +
  geom_point(data = mean_data, aes(x = mean_satisfy_index, y = vote_cat, color = vote_cat), 
             size = 4, shape = 18) +
  scale_fill_discrete(name = "Vote Category") +
  scale_color_discrete(name = "Vote Category") +
  theme_minimal() +
  labs(title = "",
       x = "BTA Satisfaction Index",
       y = "") +
  theme(legend.position = "none")
p
# ggsave("voter_bta_plot.jpg", width = 6, height = 4)

dat$maguindanao_cat <- ifelse(dat$maguindanao, "Maguindanao", "Other")
dat$maguindanao_cat <- factor(dat$maguindanao_cat, levels = c("Other", "Maguindanao"))
mean_data_maguindanao <- dat %>%
  group_by(maguindanao_cat) %>%
  summarize(mean_satisfy_index = mean(satisfy_index, na.rm = TRUE))

## Figure 11
ggplot(dat, aes(x = satisfy_index, y = factor(maguindanao_cat), fill = factor(maguindanao_cat))) +
  geom_density_ridges(scale = 0.9, alpha = 0.8) +
  geom_point(data = mean_data_maguindanao, aes(x = mean_satisfy_index, y = factor(maguindanao_cat), color = factor(maguindanao_cat)), 
             size = 4, shape = 18) +
  scale_fill_brewer(name = "Maguindanao", palette = "Pastel2") +
  scale_color_brewer(name = "Maguindanao", palette = "Set2") +
  theme_minimal() +
  labs(title = "",
       x = "BTA Satisfaction Index",
       y = "") +
  theme(legend.position = "none")
# ggsave("Maguindanao_bta_plot.pdf", width = 6, height = 4)

## Figure 12
m3q1_full <- lm(-scale(m3q1) ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                  overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                  displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(m3q1_full, exclude = exc_m3, term = fullterm, order = fullord)
# ggsave("m3q1_fullcoef.pdf", width = 6, height = 4)

## Figure 13
m3q2_full <- lm(-scale(m3q2) ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                  overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                  displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(m3q2_full, exclude = exc_m3, term = fullterm, order = fullord)
# ggsave("m3q2_fullcoef.pdf", width = 6, height = 4)

## Figure 14
m3q3_full <- lm(-scale(m3q3) ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                  overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                  displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(m3q3_full, exclude = exc_m3, term = fullterm, order = fullord)
# ggsave("m3q3_fullcoef.pdf", width = 6, height = 4)

## Figure 15
m3q4_full <- lm(-scale(m3q4) ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                  overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                  displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(m3q4_full, exclude = exc_m3, term = fullterm, order = fullord)
# ggsave("m3q4_fullcoef.pdf", width = 6, height = 4)

## Indicator of support for peace agreement (not included in the appendix)
m2q_full <- lm(-scale(m2q) ~ m2g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                 overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                 displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
exc_m2 <- c("m2gT1", "m2gT2", "m2gT3", "(Intercept)")
coefplot(m2q_full, exclude = exc_m2, term = fullterm, order = fullord)
# ggsave("m2q_fullcoef.pdf", width = 6, height = 4)

## Support for peace agreement clauses

dat$pa_sup <- dat$m1q4 + dat$m1q5 + dat$m1q6
dat$pa_sup_index <- -scale(dat$pa_sup)

## Figure 8
paindx_full <- lm(pa_sup_index ~ male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                    overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                    displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(paindx_full, exclude = "(Intercept)", term = fullterm, order = fullord)
# ggsave("paindx_fullcoef.jpg", width = 6, height = 4)



exc_m1 <- c("m1q4gT1", "m1q4gT2", "m1q4gT3", "m1q4gT4", "m1q4gT5", "m1q4gT6", "m1q5gT1", "m1q5gT2", "m1q5gT3",
            "m1q5gT4", "m1q5gT5", "m1q5gT6", "m1q6gT1", "m1q6gT2", "m1q6gT3", "m1q6gT4", "m1q6gT5", "m1q6gT6",
            "(Intercept)")

## Figure 17
pa1_full <- lm(-scale(m1q4) ~ m1q4g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                 overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                 displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(pa1_full, exclude = exc_m1, term = fullterm, order = fullord)
# ggsave("pa1_fullcoef.pdf", width = 6, height = 4)

## Figure 18
pa2_full <- lm(-scale(m1q5) ~ m1q5g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                 overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                 displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(pa2_full, exclude = exc_m1, term = fullterm, order = fullord)
# ggsave("pa2_fullcoef.pdf", width = 6, height = 4)

## Figure 19
pa3_full <- lm(-scale(m1q6) ~ m1q6g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                 overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                 displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
coefplot(pa3_full, exclude = exc_m1, term = fullterm, order = fullord)
# ggsave("pa3_fullcoef.pdf", width = 6, height = 4)

### Regression on priority political issue in Bangsamoro
## Redefine the issues
dat$issue_ml <- dat$pol_issue
## Less popular issues are aggregated into "other"
dat$issue_ml[dat$issue_ml %in% c("disaster", "health", "land", "education")] <- "other"
dat$issue_ml <- factor(dat$issue_ml, levels = c("other", "economy", "peace"))


issue_reg <- multinom(issue_ml ~ male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                        overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant + peaceful + witness_reb + witness_rido +
                        displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
stargazer(issue_reg)
# Extract coefficients
coefficients <- coef(issue_reg)

# Calculate confidence intervals
conf_intervals <- confint(issue_reg)

# Transform coefficients and confidence intervals to oddat ratios
odds_ratios <- exp(coefficients)
conf_intervals_exp <- exp(conf_intervals)

# Prepare a data frame for plotting
terms <- rep(colnames(coefficients), each = nrow(coefficients))
responses <- rep(rownames(coefficients), times = ncol(coefficients))

coef_df <- data.frame(
  term = terms,
  response = responses,
  estimate = as.vector(coefficients),
  conf.low = as.vector(t(conf_intervals[, 1,])),
  conf.high = as.vector(t(conf_intervals[, 2,]))
)

coef_df <- coef_df[coef_df$term != "(Intercept)",]

coef_df$term <- rep(fullterm, each = 2)
coef_df$term <- factor(coef_df$term, levels = rev(fullterm))

# Plot odds ratios with confidence intervals
## Figure 9
ggplot(coef_df, aes(x = term, y = estimate, color = response)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(width = 0.5)) +
  scale_color_brewer(palette = "Dark2") +
  coord_flip() +
  theme_minimal() +
  labs(title = "",
       x = "",
       y = "Coefficients",
       color = "Response Category") +
  geom_hline(yintercept = 0, linetype = "dashed")
# ggsave("mlogit.jpg", width = 6, height = 6)


########## Immigrartion, further analysis

indx_im1_full <- lm(satisfy_index ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                      overhighschool + maguindanao + tagalog + iranun + urban + cotcity + immigrant_5 + immigrant_10 +
                      immigrant_long + peaceful + witness_reb + witness_rido + displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
summary(indx_im1_full)

indx_im2_full <- lm(satisfy_index ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                      overhighschool + maguindanao + tagalog + iranun + urban + cotcity +
                      othermun + otherprov + noBARMM + peaceful + witness_reb + witness_rido +
                      displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
summary(indx_im2_full)

indx_im3_full <- lm(satisfy_index ~ m3g + male + richer_median_pc + old + homemaker + selfemploy + farmer + worker + overelement +
                      overhighschool + maguindanao + tagalog + iranun + urban + cotcity +
                      (othermun + otherprov + noBARMM) + (immigrant_5 + immigrant_10 + immigrant_long) +
                      peaceful + witness_reb + witness_rido + displace + clan_affiliation + vote + UBJP_vote + better_knowledge, dat)
summary(indx_im3_full)

# Extract coefficients and CIs from the models
coef_im1 <- tidy(indx_im1_full, conf.int = TRUE)
coef_im2 <- tidy(indx_im2_full, conf.int = TRUE)

# Filter for the desired variables
vars1 <- coef_im1 %>% filter(term %in% c("immigrant_5", "immigrant_10", "immigrant_long"))
vars2 <- coef_im2 %>% filter(term %in% c("othermun", "otherprov", "noBARMM"))

# Add a model identifier to each data frame
vars1$model <- "Temporal"
vars2$model <- "Spatial"

# Combine the data frames
combined_vars <- rbind(vars1, vars2)

# Create a custom x-axis
combined_vars$custom_x <- factor(combined_vars$term, levels = c("immigrant_5", "immigrant_10", "immigrant_long", "othermun", "otherprov", "noBARMM"))

# Assign numeric values to custom_x for plotting
combined_vars$custom_x_num <- as.numeric(combined_vars$custom_x)

# Figure 7
ggplot(combined_vars, aes(x = custom_x_num, y = estimate, ymin = conf.low, ymax = conf.high, color = model)) +
  geom_pointrange(position = position_dodge(width = 1), na.rm = TRUE) +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6), labels = c("~ 5 years", "5-10 years", "10 ~ years", "In Maguindanao", "In BARMM", "Outside BARMM")) +
  labs(title = "",
       x = "Variables",
       y = "Coefficients",
       color = "Model") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  geom_hline(yintercept = 0, lty = 2) +
  scale_color_brewer(palette = "Dark2")
# ggsave("immigrant.jpg", width = 6, height = 4)


####### Clan affiliation, additional analysis

clan_full <- lm(clan_affiliation ~ male + richer_median_pc + old + homemaker + selfemploy + farmer + worker +
                  overelement + overhighschool + maguindanao + tagalog + iranun +
                  urban + cotcity + peaceful + witness_reb + witness_rido + immigrant +
                  displace + vote + UBJP_vote + better_knowledge, dat)

## Figure 10
coefplot(clan_full, exclude = exc_m3, term = c("Male", "Richer", "Older", "Homemaker", "Selfemploy",
                                               "Farmer", "Worker", "ElementarySchool", "Highschool",
                                               "Maguindanao", "Tagalog", "Iranun", "Urban", "Cotabato City",
                                               "Immigrant", "Peaceful", "Witness Insurgency", "Witness Rido",
                                               "Displacement", "Voted", "Voted for UBJP", "Political Knowledge"))
# ggsave("clan_reg.pdf", width = 6, height = 4)
